Journal of Computational Chemistry

RNA function is intrinsically linked to its structural polymorphism, with molecules exploring the heterogeneous conformational ensembles resulting from complex energy landscapes. These landscapes arise from competing interactions, small energetic separations between microstates, and strong coupling to the environment, posing significant challenges for both experimental characterization and molecular simulation. In this chapter, we review current computational strategies that aim to explore RNA conformational ensembles in silico, with a specific focus on energy landscape-based approaches and atomistic simulations. We discuss key limitations related to sampling efficiency, force-field accuracy, and ensemble analysis, and illustrate their impact through case studies on a self-cleaving ribozyme and an H-type pseudoknot. Finally, we highlight emerging directions, including closer integration with experimental data and the growing role of machine learning, which will probably reinforce the predictive power of in silico RNA energy landscape exploration.

Transition metal based compounds are promising therapeutic agents, particularly in cancer treatment. However, predicting their binding sites remains a major challenge. In this work, we investigate the applicability of two tools, Metal3D and Metal1D, for this purpose. Although originally trained to predict zinc ion binding sites only, both predictors successfully identify several experimentally observed binding sites for transition metal complexes directly from apo protein structures. At the same time, we highlight current limitations, such as the sensitivity to side-chain conformations, and discuss possible strategies for improvement. This work provides a first step toward establishing a robust computational pipeline in which rapid and low-cost predictors are able to identify putative hotspots for transition metal binding, which can then be refined using more accurate but computationally demanding methods. Author summaryTransition metals play a crucial role as therapeutic agents, especially in cancer therapy. However, the prediction of their binding site locations is challenging, as accurate computational methods often require time-consuming simulations, making them impractical when many possible binding sites must be explored. In this work, we explored the capability of two binding site predictors, originally developed to locate metal ions in proteins, to identify binding sites for more complex covalently-bound transition metal based agents. We found that these tools can often identify the experimentally-known binding regions, even when starting from the apo structure, in which the protein does not already contain the metal compound. At the same time, our results show clear limitations in more challenging cases, particularly when the binding involves only a single amino acid or when the binding site undergoes major structural rearrangements upon binding. Overall, our study shows that fast predictors can provide valuable early insights in the investigation of the binding sites of covalently-bound transition metal based compounds. When combined with more accurate simulation techniques, they can help focus computational efforts and ultimately support the rational design of transition metal based drugs.

This paper describes MMP, a three-stage framework for systematic quantum optimization of constrained molecular docking problems. The protocol addresses the "formulation bottleneck"--the critical challenge of translating constrained optimization problems into valid QUBO (Quadratic Unconstrained Binary Optimization) formulations for quantum solvers. MMP replaces heuristic penalty tuning with data-driven calibration through: (1) classical solution-space analysis to validate fragment libraries before quantum deployment, (2) systematic penalty sweeps to identify optimal "Goldilocks Zone" coefficients, and (3) MAC-QAOA (MMP Adaptive Constraint QAOA) with layer-dependent penalty decay. Preliminary benchmarks on synthetic constrained optimization problems demonstrate 99.7% solution validity at identified elbow points and 25.5% improvement in solution quality over static-penalty QAOA. MMP is hardware-agnostic but designed for near-term devices including Pasqals Orion Gamma (140+ qubits). The theoretical framework, algorithmic details, and preliminary validation results of the protocol are discussed, establishing a systematic methodology for quantum-augmented optimization workflows for drug discovery. All benchmarks are conducted on synthetic constrained optimization instances that reproduce structural features of docking formulations; application to real molecular docking targets is left for future work.

Recently, a number of tools have been released that generate ensembles of protein structures based on artificial intelligence (AI) approaches. Although ensembles generated by the tools differ significantly, we demonstrate a computational path to harmonizing the various outputs under a stationary condition using two complementary physics-based approaches. In the first stage, the AI ensemble is used to seed a weighted ensemble (WE) simulation, promoting relaxation toward the steady state. In the second stage, trajectory segments generated by WE are reweighted to steady state using the recently developed RiteWeight (RW) algorithm. We applied this approach to generate an atomically-detailed equilibrium ensemble of unliganded adenylate kinase conformations, starting from ensembles produced by three AI tools: AFSample2, ESMFlow-PDB (trained from PDB structures), and ESMFlow-MD (trained from molecular dynamics simulation data). Dramatic differences in the AI-generated ensembles are largely erased during the WE-RW process, yielding a consistent description of the equilibrium ensemble for a given force field.

The incorporation of non-canonical amino acids (ncAAs) is a powerful strategy for introducing novel chemical functions into proteins. Molecular dynamics (MD) simulations are essential for understanding the structural and dynamic effects of these modifications, yet the creation of accurate force field parameters for ncAAs remains a significant bottleneck. Current parameterisation methods are often inaccurate or computationally expensive. To address this, we present drFrankenstein, an automated pipeline for generating AMBER force field parameters for ncAAs. drFrankenstein is a robust and accessible tool that streamlines the parameterisation workflow, enabling the routine use of MD simulations to study the behaviour of ncAA-containing proteins.

An unresolved challenge in the field of computational protein design is to create proteins that bind non-protein partners, e.g. DNA, RNA, and small molecules. Most machine learning (ML) algorithms for protein design can only work with systems composed entirely of amino acids, and therefore cannot be directly applied to this task. The few algorithms that accommodate non-proteins still represent amino acids differently than other molecules, and therefore cannot easily recognize the similarity between a sidechain and a small molecule that share a functional group. We introduce a new method, called AtomPaint, that avoids these limitations by employing a fully-atomic representation of protein structure. Starting from a model of a desired binding interaction, our method proceeds by (i) converting that model to a 3D image, (ii) masking out the parts of that image that need to be redesigned, (iii) using a diffusion model to inpaint the masked voxels, then (iv) using a classification model to identify the amino acids in the inpainted image. Both models are SE(3)-equivariant ResNets, and were trained on a dataset of structures from the Protein Data Bank (PDB) curated to emphasize protein/non-protein interactions. In a sequence recovery benchmark, AtomPaint performed better than random guessing, suggesting that it understands some aspects of molecular structure. We discuss possible avenues of improvement, in the hopes that the advantages of our novel image-based approach can be fully realized.

There is a growing demand for molecular dynamics simulations to explore longer timescale behavior of giant protein-DNA complexes such as chromatin. To address this need, we extended OpenCafeMol, a GPU-accelerated residue-level coarse-grained molecular dynamics simulator originally developed for proteins and lipids, to support 3SPN.2 and 3SPN.2C DNA models. We also implemented a hydrogen-bond-type many-body potential to model DNA-protein interactions more accurately. To further improve computational efficiency, we introduced a localized scheme for calculating base-pairing and cross-stacking interactions. Benchmark tests show that OpenCafeMol on a single GPU achieves up to 200-fold speed-up for DNA-only systems and up to 100-fold speed-up for DNA-protein complexes compared to CPU-based simulations. To demonstrate the capability of our implementation for long-timescale biological processes, we simulated an archaeal SMC-ScpA complex undergoing DNA translocation via segment capture (a proposed mechanism for DNA loop extrusion) in the presence of a DNA-bound obstacle. We observed continuous captured-loop growth accompanied by obstacle bypass within the segment capture framework.

Millions of experimental and AI-predicted protein structures are now available, and the biosynthetic promise of bespoke proteins is increasingly within reach. The functional characterization challenge thus posed cannot be addressed by experimental techniques alone. Molecular dynamics (MD) simulations offer functional screening with atomic resolution, yet accessibility remains limited. Existing computational chemistry software presents stark trade-offs whereby powerful tools require extensive expertise and manual effort, or user-friendly programs function as black boxes that obscure critical preparation decisions. Herein, we present ProPrep, an interactive workflow manager that guides users through expert-quality MD preparation by showing the what, why, and how of each step while automating tedious manual operations. Within a single workspace, ProPrep integrates (1) downloading structures from multiple sources (PDB, AlphaFold, AlphaFill), (2) performing homology searches, (3) aligning structures, (4) curating and repairing structural issues, (5) applying mutations, (6) parameterizing specialized residues, (7) converting redox-active sites to forcefield-compatible forms, (8) generating topology and coordinate files, and (9) configuring, executing, and analyzing simulations with active monitoring of key quantities via ASCII visualizations. A key innovation is ProPreps extensible transformer framework for detecting, defining, and transforming redox-active sites--including mono- and polynuclear metal centers, organic cofactors, and redox-active amino acids--for forcefield compatibility. We demonstrate the full workflow on a 64-heme cytochrome nanowire bundle (PDB: 9YUQ), proceeding from a PDF file to energy minimization of the solvated system (467,635 atoms) for constant pH molecular dynamics--a process demanding 4,819 PDB record modifications and 610 bond definitions--in 18 minutes of user interaction. The entire process is recorded in an interactive session log that can be shared and replayed for reproducibility, making simulation setup a fully transparent process that relies on what was done instead of what was remembered and reported.

As high-performance computing provides the ability to generate and analyze ever larger simulation trajectories, the challenges in learning, applying, and sharing the best analytical practices become more salient. Extracting reproducible scientific insights from simulation requires a thorough understanding of many computing topics unrelated to the molecular systems being modeled and simulated. While the rapid development of the technologies used for analysis makes previously impossible studies into routine work, the growing repertoire of software combined with the specificity of the ecosystems that they rely on can easily break the programs used in older studies. In this work, we present ST-Analyzer, a simulation trajectory analysis suite with command-line (CLI) and graphical (GUI) user interfaces. ST-Analyzer is distributed freely as an open-source conda-forge package with support for macOS, Linux, and Windows (via WSL2). Besides facilitating several common analysis tasks, the GUI shows users the exact commands necessary to repeat the same tasks on the command-line. We demonstrate ST-Analyzers capabilities by reproducing several results from previously published simulation studies on the lipid parameters of heterogeneous biomembranes and the behavior of a SARS-CoV-2 spike protein-antibody complex. We expect ST-Analyzer to be useful to experts for quickly setting up common analysis tasks and to nonexperts as a guided introduction to simulation analysis using both GUI and CLI. ST-Analyzer is freely available at https://github.com/nk53/stanalyzer.

The unfolding or melting temperature (TM) is a central quantity to characterize the stability of proteins and other biopolymers. The accurate prediction of protein melting temperatures by molecular mechanics force field simulations is highly desirable for many biophysical and biotechnological applications. Since the time scales for protein (un-)folding are hardly accessible in conventional MD (cMD) simulations, enhanced sampling techniques such as Temperature Replica Exchange Molecular Dynamics (TREMD) are typically employed. However, TREMD simulations are computationally very demanding especially if large temperature ranges need to be covered. Additionally, if the TM is initially unknown, setting up TREMD simulations is often challenging. To find the optimal initial conditions for such simulations, we describe their performance based on a theoretical model, which we validate on a minimalistic Markov Chain Monte Carlo (MCMC) simulation setup. In an effort to reduce the computational demand, we have investigated the possibility to use small sets of TREMD temperature ladders placed iteratively in the vicinity of a TM estimate. Different TREMD setups were extensively tested on the fast-folding protein Chignolin. We found that appropriate starting conformations lead to significantly faster convergence. Furthermore, we found that, in practice, combining multiple small temperature ladders can be advantageous in comparison to one single temperature ladder. Based on our findings, we formulate practical recommendations on how to set up TREMD for protein melting with optimal efficiency.

Molecular dynamics simulation is a powerful computational technique used for predicting and understanding the dynamic behavior of biomolecular systems. Steered molecular dynamics (SMD) simulations enable the study of force-induced processes in biomolecules, effectively mimicking single-molecule force spectroscopy experiments probing protein unfolding and receptor-ligand unbinding. Given the stochastic nature of these mechanical events, accurately exploring the dynamic behavior of biomolecules and extracting accurate physical information requires several in-silico experiments. This includes performing many pulling simulations at different velocities or force loading rates. The large amount of data obtained from these simulation sets requires efficient automated data processing tools. We present PySteMoDA, a novel Python package with a user-friendly graphical interface specifically designed for constant-velocity SMD data analysis. The automated force peak detection methods reduce user bias, improve accuracy, and accelerate data analysis. The package also allows identification of residues involved in mechanical events through computation of the time-dependent mechanical work and correlation factors between residue pairs. This package not only addresses automated data processing in SMD simulations and accurate parameter extraction, but also significantly enhances accessibility and usability. Through PySteMoDA, users can efficiently analyze simulation data without the barrier of coding, facilitating a wider range of investigations and insights in the field of computational biochemistry and biophysics.

Accessing deeply buried binding sites remains a major challenge in structure-based drug discovery, where accurate description of both protein dynamics and ligand binding pathways is required. Funnel metadynamics enables simulation of complete binding processes but is computationally demanding at the all-atom resolution. By adopting the Martini 3 force field, coarse-grained funnel metadynamics (CG-FMD) substantially reduces computational requirements while retaining enhanced sampling capabilities. In this work, we assess the capability of CG-FMD to model ligand recognition at the deeply buried colchicinoids site of the tubulin {beta}-heterodimer, a multisite protein of strategic importance. We investigated the binding of colchicine, podophyllotoxin and combretastatin-A4, recovering free energy profiles with improved statistical convergence compared to AA-FMD and comparable to experimental references. In particular CG-FMD binding free energies present mean absolute errors between 3 and 10 kJ mol-1. These results propose CG-FMD as an efficient, physics-based framework for probing ligand binding to challenging sites.

AO_SCPLOWBSTRACTC_SCPLOWThe Martini coarse-grained force field is widely used for biomolecular simulations by a large and rapidly expanding community worldwide. Over time, the development of Martini parameters, tools, and documentation has become increasingly dispersed across numerous research groups, leading to fragmentation and making it challenging for users and developers to keep track of the latest models, software, and best practices. Consequently, the development of Martini as a genuinely community-driven process has grown into a bottleneck. In response, the Martini Force Field Initiative (MFFI) has been established as an open-science effort to coordinate and support the collaborative development of all Martini resources. Here, we introduce the MFFI web portal, a platform designed around five core pillars: (i) avoiding reliance on a single group or local server; (ii) minimizing long-term maintenance overhead; (iii) reducing technical barriers for contributions; (iv) providing a unified home for parameters, tools, tutorials, example workflows, and research outputs; and (v) enabling timely dissemination of updates to the community. To achieve this, we use Quarto to generate a static website authored in Markdown, lowering the technical barrier to making contributions, and serverless architectures on Amazon Web Services for scalable, event-triggered backend operations. The source code is hosted in a public GitHub repository under an MIT license, with automated deployment via GitHub Actions and a contribution model based on pull requests for quality control. This design creates a sustainable, low-maintenance, and collaborative infrastructure that consolidates Martini resources and supports transparency. More broadly, our design exemplifies a transferable pattern for building open, community-oriented platforms for molecular modeling and computational science.

The MArtini Database (MAD) web server (https://mad.ens-lyon.fr/) has been updated with new tools, models, and capabilities to support a broader range of molecular systems for the Martini coarse-grained (CG) force field. The most notable addition is the MAD:Polymer Builder, enabling the automated construction of CG polymer models with varying architectures and complexities. The server now incorporates the latest developments in Martini 3, providing enhanced control within the MAD:Molecule Builder over the conversion of all-atom structures into CG representations, including the implementation of G[o]Martini 3 for protein complexes and water-protein interaction biases. The MAD:Polymer Builder and MAD:Molecule Builder are both adapted to work with intrinsically disordered proteins and domains. Substantial progress has been made in expanding the MAD:Database, making a growing library of Martini-ready compounds readily accessible across the entire MAD ecosystem. These advances position MAD as a comprehensive and evolving platform for the preparation of diverse systems in CG molecular simulations.

With the advance of hardware and software for molecular dynamics simulation it has become routine to obtain trajectories that are tens of microseconds in duration for all kinds of protein machinery. This shifts the burden of work onto analysis of the simulation data and opens opportunities for more rigorous and reproducible observations on mechanism. Toward this end we developed an investigator-blind analysis pipeline which operates on featurized simulation data, performs unsupervised clustering, and then identifies which input features are most discriminatory of cluster identity. Application of this pipeline to a large set of G-protein coupled receptor simulation data shows that it identifies several well-known microswitches. Inspection of these structural elements reveals changes in conformation that are known to accompany functional transitions of the receptor. In addition to these known structural elements the analysis also identifies two possibly new structural motifs: the kink in transmembrane helix 2, and a coupled "piston-like" motion of TM2 and TM3.

We present DIVIsive N-ary Ensembles (DIVINE), a deterministic, top-down clustering framework designed for molecular dynamics (MD) trajectories. DIVINE constructs a complete clustering hierarchy by recursively splitting clusters based on n-ary similarity principles, avoiding the need for O(N2) pairwise distance matrices. It supports multiple cluster selection criteria, including a weighted variance metric, and deterministic anchor initialization strategies such as NANI (N-ary Natural Initiation), ensuring reproducible and structurally meaningful partitions. Testing DIVINE up to a 305 s folding trajectory of the villin headpiece (HP35) revealed that it matched or exceeded the clustering quality of bisecting k-means while reducing runtime and eliminating stochastic variability. Its single-pass design enables efficient exploration of clustering resolutions without repeated executions. By combining scalability, interpretability, and determinism, DIVINE offers a robust and practical alternative to conventional MD clustering methods. DIVINE is publicly available as part of the MDANCE package: https://github.com/mqcomplab/MDANCE.

Molecular dynamics (MD) simulations are a powerful tool for investigating biomolecular dynamics underlying biological functions. However, the accessible spatiotemporal scales of conventional all-atom simulations remain limited by high computational costs. Coarse-graining reduces these costs by decreasing the number of interaction sites and enabling longer timesteps. In extreme cases, proteins are represented as single spherical particles; while such approximations facilitate cellular-scale simulations, they often sacrifice essential structural information, such as molecular shape and interaction anisotropy. Here, we present CGRig, a rigid-body protein model with residue-level interaction sites designed for long-time, large-scale simulations. In CGRig, each protein is treated as a single rigid-body embedding residue-level interaction sites. Its translational and rotational motions are described by the overdamped Langevin equation incorporating a shape-dependent friction matrix. Intermolecular interactions are calculated using G[o]-like native contact potentials, Debye-Huckel electrostatics, and volume exclusion. We validated that CGRig accurately reproduces the translational and rotational diffusion coefficients expected from the friction matrix for an isolated protein. For dimeric systems, the model successfully maintained native complex structures. Furthermore, two initially separated proteins converged into the correct complex with an association rate consistent with all-atom simulations. Notably, CGRig achieved a simulation performance exceeding 17 s/day for a 1,024-molecule system. These results demonstrate that CGRig provides an efficient framework for simulating protein assembly while retaining residue-level interaction specificity, making it a valuable tool for investigating large-scale biomolecular self-assembly.

Kinase family proteins constitute the second largest protein class targeted in drug development efforts, most prominently to treat cancer, but also several other diseases associated with kinase dysfunction. In this work we focus on type II kinase inhibitors which bind to the "classical" inactive conformation of the protein kinase catalytic domain where the DFG motif has a "DFG-out" orientation and the activation loop is folded. Many Tyrosine kinases (TKs) exhibit strong binding affinity with a wide spectrum of type II inhibitors while serine/threonine kinases (STKs) often bind more weakly. Recent work suggests this difference is largely due to differences in the folded to extended conformational equilibrium of the activation loop between TKs vs. STKs. The binding affinity of a type II inhibitor to its kinase target can be decomposed into a sum of two contributions: (1) the free energy cost to reorganize the protein from the active to inactive state, and (2) the binding affinity of the type II inhibitor to the inactive kinase conformation. In previous work we used a Potts statistical energy potential based on sequence co-variation to thread sequences over ensembles of active and inactive kinase structures. The threading function was used to estimate the free energy cost to reorganize kinases from the active to classical inactive conformation, and we showed that this estimator is consistent with the results of molecular dynamics free energy simulations for a small set of STKs and TKs. In the current study, we analyze the results of a large-scale study of the binding affinities of 50 type II inhibitors to 348 kinases, of which the results for 16 of the 50 type II inhibitors were reported in an earlier study (the "Davis dataset"); the binding data for the remaining 34 type II inhibitors to the panel of 348 kinases were recently obtained (the "Schrodinger dataset"). We use the Potts statistical energy model to investigate the contribution of protein reorganization to the selectivity of the large kinase panel against the set of 50 type II inhibitors, and find that protein reorganization makes a significant contribution to the selectivity. The AUC of the receiver-operator characteristic curve is [~]0.8. We report the results of an internal "blind test", that shows how Potts threading energies can provide more accurate estimates of kinase selectivity than corresponding predictions using experimental results of small sample size. We discuss why two STK phylogenetic kinase families, STE and CMGC, appear to contain many outliers, and how to improve the ability to predict kinase selectivity with a more complete analysis of the kinase conformational landscape. We compare the performance of Potts threading for predicting binding properties of the large set of (50) Type II inhibitors to 348 kinases, with those of a sequence-based purely machine learning model, DeepDTAGen, a publicly available machine learning model that was trained on the complete Davis dataset, including both Type I and Type II kinase inhibitors. We observe that DeepDTAGen performs well on binding predictions for the 16 type II inhibitors in the Davis dataset, but performs poorly on binding predictions for the 34 type II inhibitors against 348 kinases in the Schrodinger dataset.

Phase separated condensates are recognized as a ubiquitous mechanism of spatial organization in cell biology. Biophysical modeling of condensates provides critical insights into the dynamics and functions of these subcellular structures that are difficult to extract via experiments. Here we present an efficient computational pipeline, CASPULE (Condensate Analysis of Sticker Spacer Polymers Using the LAMMPS Engine), to simulate and analyze the biological condensates made of sticker-spacer polymers. CASPULE implements a unique force field that combines traditional Langevin dynamics with a "detailed balance proof" protocol for single-valent bond formation between stickers. This framework allows us to study the non-trivial biophysics that emerge out of the single-valent sticker interactions coupled with the effect of separation in energetic contribution by stickers and spacers. We provide detailed documentation on how to setup the simulation environment, perform simulations and analyze the results. Through case studies, we highlight the utility and efficacy of our pipeline. Importantly, we provide statistical parameters to characterize the cluster size distribution often observed in biological systems. We envision this tool to be broadly useful in decoding the interplay of kinetics and thermodynamics underlying the formation and function of biological condensates.

This study presents a comprehensive analysis of the dynamic properties and allosteric regulation mechanisms of Class A G protein-coupled receptors (GPCRs) by integrating molecular dynamics (MD) simulations with nuclear magnetic resonance (NMR) relaxation measurements. Utilizing generalized order parameters derived from NMR data and MD trajectories, we quantitatively assess conformational entropy changes that occur during receptor activation and ligand binding events. This approach enables a detailed characterization of protein flexibility at multiple timescales, revealing how dynamic fluctuations contribute to allosteric signal transmission within the receptor. Our results demonstrate that conformational entropy plays a pivotal role in modulating the functional states of Class A GPCRs, influencing the equilibrium between inactive and active conformations. By elucidating the interplay between structural dynamics and allostery, this work advances the molecular-level understanding of GPCR function and highlights the importance of entropy-driven effects in receptor signaling. The integrative methodology and findings provide a valuable framework for future investigations aimed at targeting receptor dynamics in drug discovery and rational design of allosteric modulators.